!> 预估校正法, 用于求解非线性系统
module seakeeping_predictor_corrector

    use seakeeping_kinds, only: wp => sk_real_kind
    implicit none

    private
    public :: predictor_corrector

    abstract interface
        !> `fcn` 计算 ODE 的导数
        subroutine fcn(t, y, yp)
            import wp
            real(wp), intent(in) :: t       !! 自变量
            real(wp), intent(in) :: y(:)    !! 因变量
            real(wp), intent(out) :: yp(:)  !! 导数
        end subroutine fcn
    end interface

contains

    !> 预估校正法
    subroutine predictor_corrector(func, t, y, dt, m, n)
        procedure(fcn) :: func                                  !! 回调函数
        real(kind=wp), intent(inout) :: t                       !! 时间
        integer, intent(in) :: n                                !! 长度
        real(kind=wp), intent(inout), dimension(*) :: y         !! 因变量
        real(kind=wp), intent(in) :: dt                         !! 时间步长
        integer, intent(in) :: m                                !! 积分次数
        real(wp), dimension(n) :: y_pred, k1, k2
        real(wp) :: h2
        integer :: i

        h2 = dt/2
        do i = 1, m
            call func(t, y(:n), k1)         ! 预测
            y_pred(:n) = y(:n) + k1*dt
            t = t + dt
            call func(t, y_pred(:n), k2)    ! 校正
            y(:n) = y(:n) + (k1 + k2)*h2
        end do

    end subroutine predictor_corrector

end module seakeeping_predictor_corrector
